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Critical point of Nf = 3 QCD from lattice simulations in the canonical ensemble 
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A canonical ensemble algorithm is employed to study the phase diagram of Nf — 3 QCD using 
lattice simulations. We lock in the desired quark number sector using an exact Fourier transform of 
the fermion determinant. We scan the phase space below T c and look for an S-shape structure in 
the chemical potential, which signals the coexistence phase of a first order phase transition in finite 
volume. Applying Maxwell construction, we determine the boundaries of the coexistence phase at 
three temperatures and extrapolate them to locate the critical point. Using an improved gauge 
action and improved Wilson fermions on lattices with a spatial extent of 1.8 fm and quark masses 
close to that of the strange, we find the critical point at Te = 0.925(5) T c and baryon chemical 
potential = 2.60(8) T c . 
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QCD is expected to have a rich phase diagram at finite 
temperature and finite density. Current lattice calcula- 
tions have shown that the transition from the hadronic 
phase to QGP phase is a rapid crossover [HE]- For large 
baryon chemical potential and very low temperature, a 
number of models suggest that the transition is a first 
order. If this is the case, when the chemical potential 
is lowered and temperature raised, this first order phase 
transition is expected to end as a second order phase tran- 
sition point — the critical point. However, lattice QCD 
simulations with chemical potential are difficult due to 
the notorious "sign problem" . The majority of current 
simulations are focusing on small chemical potential re- 
gion n q /T -C 1 where the "sign problem" appears to be 
under control. Up to now, all the Nf = 3 or 2+1 sim- 
ulations are based on the grand canonical ensemble (T, 
fiB as parameters) with staggered fermions. The results 
from the multi-parameter reweighting [3], Taylor expan- 
sion with small n [H [5] and the curvature of the critical 
surface [Bj are not settled and need to be cross-checked. 
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FIG. 1. Conjectured QCD phase diagram 
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Even the existence of the critical point is in question [BJ . 
We employ an algorithm, which is not restricted to small 
chemical potential because of the mitigation of the sign 
problem under the current parameter settings, to study 
this problem. 

In this letter, we adopt an exact Monte Carlo algo- 
rithm [THS] based on the canonical partition function [TOl - 
IT5] which is designed to alleviate the determinant fluctu- 
ation problems. As it turns out, the sign fluctuations are 
not serious on the lattices used in the present study, as we 
shall see later. In the canonical ensemble simulations in 
finite volume, the coexistence phase of a first order phase 
transition has a characteristic S-shape as a function of 
density due to the surface tension. This finite-volume 
property has been exploited successfully to identify the 
phase boundaries via the Maxwell construction in studies 
of phase transition with the staggered fermions [HI [TS] 
and clover fermions [TBJ for the Nf — 4 case which is 
known to have a first order phase transition at /i = 0. In 
these benchmark studies the boundaries were identified 
at three temperatures below T c , and they were extrap- 
olated in density and temperature to show that the in- 
tersecting point indeed coincides with the independently 
identified first order transition point at T c and /x = [TBJ . 
In view of the success of the Nf = 4 study, we extend this 
method to the more realistic Nf — 3 case [17j [18] . Al- 
though the real world contains two light quarks and one 
heavier strange quark, the three degenerate flavor case 
has a similar phase structure. Our primary goal in this 
study is to determine whether a first order phase tran- 
sition exists for Nt — 3 and where the critical point is 
located. 

With the aid of recently developed matrix reduction 
technique [191121) . we scan the chemical potential as a 
function of baryon number for four temperatures below 
T c which is determined at zero chemical potential, and we 
observe clear signals for a first order phase transition for 
temperatures below 0.93 T c . The phase boundaries of the 
coexistence phase are determined and then extrapolated 
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in temperature and density to locate the critical point 
at T E = 0.927(5) T c and /if = 2.60(8) T c . Our results 
are based on simulations on 6 3 x 4 lattices with clover 
fermion action with quark masses which correspond to 
the pion mass from 750 MeV for the lowest temperature 
to 775 MeV for the highest temperature. 

The canonical partition function in lattice QCD can be 
derived from the fugacity expansion of the grand canon- 
ical partition function, 



Z(V,T,ri = Y i Z c (V,T,k)e 



Hk/T 



k 
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where k is the net number of quarks (number of quarks 
minus the number of anti-quarks) and Zq is the canonical 
partition function. Using the fugacity expansion, it can 
be shown that the canonical partition function can be 
written as a Fourier transform of the grand canonical 
partition function, 
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upon introducing an imaginary chemical potential /j, = 
i(j>T. After integrating out the fermionic part in Eq. ([2]), 
we get an expression 



Z c (y,T,k)= V\Je- s <> iu Uet k M N f(U), (3) 



where 



det k M N f (U) 
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d<£e- ifc * det M(m,ct);U) N f, 



(4) 

is the projected determinant for the fixed net quark num- 
ber k. Nf is the number of flavors. We shall use the re- 
cently developed matrix reduction technique to compute 
the projected determinant exactly [T§] . 

Using charge conjugation symmetry, one can show that 
defy-M f(U) is real, but not necessarily positive. Due 
to the sign fluctuation, there can potentially be a sign 
problem at large quark number and low temperature. 
For more detailed discussion about the properties of the 
canonical ensemble, we refer the reader to Ref. [IB]. To 
simulate Eq. ^ dynamically, we rewrite canonical par- 
tition function as 

Z c {V,T,k) = J V\]e- s o {u) &etM N f{U)W{U)a{U),{h) 
where 



_ |Rcdet fc M^(l/)| 
K ' detM N f(U) ' 
a(U) = Sign(Redet fe M^ ([/)). 



(6) 



Our strategy to generate an ensemble is to employ 
Metropolis accept/reject method based on the weight 
W(U) and fold the phase factor a(U) into the mea- 
surements. In short, during the simulation, the candi- 
date configuration is "proposed" by the standard Hybrid 
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FIG. 2. Average sign as a function of ub for highest and 
lowest temperatures (0.93 T c and 0.85 T c ) used in this study. 
Dashed lines represent the phase boundaries of the coexis- 
tence phase. 



Monte Carol algorithm and then an accept /reject step is 
used for the correct probability. Note the two-step sim- 
ulation with HMC and accept/reject based on W(U) re- 
duces the fluctuation problem [5] and accept /reject step 
based on the exact projected determinant detfeM^ (U) 
ensures that the simulation remains in the specific canon- 
ical sector with quark number k ^= 0. 

The lattice spacing and the pion mass are determined 
by using dynamically generated ensembles on I2 3 lattices 
for each j3. To locate the pseudo critical temperature T c , 
we varied /3 to look for the peak of the Polyakov loop 
susceptibility. We run simulations for five different vol- 
umes (6 3 , 8 3 , 10 3 , 12 3 , 16 3 x 4) and found that the peak of 
the susceptibility hardly depends on the volume. This is 
consistent with the finding on large volumes and physical 
quark masses that the finite temperature transition for 
the Nf = 3 case is a crossover at zero chemical poten- 
tial [HE]. 

To determine the location of the phase transition at 
non-zero baryon density, we pick four temperatures be- 
low T c (0.85 T c , 0.87 T c , 0.90 T c , 0.93 T c ) and vary the 
net quark number from 3 to 54 in steps of 3 (for frac- 
tional baryon number the partition function vanishes). 
This corresponds to the baryon number n B from 1 to 18 
and a density between that of the nuclear matter and 18 
times of that. The chemical potential is calculated and 
plotted as a function of the net baryon number n B . In 
the canonical ensemble, the baryon chemical potential is 
calculated by taking the difference of the free energy after 
adding one baryon, i.e. 



where 



F(n B + 1) - F{n B ) 
(n B + 1) - n B 



.iln« (7) 



Redet 3ni3+ 3M"/(^) 
|Redet3„ B M"/(?7)| ' 
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is measured in the ensemble with n B baryon number and 
() D in Eq. Q stands for the average over the ensemble 
generated with the measure |Redet3 ns M n/ (U)\. 

As a first check, we examine the magnitude of the sign 
fluctuations. The average sign in Eq. ([§ appears in the 
denominator of Eq. ^ and can lead to a sign problem 
when its error bars overlap with zero. This quantity is 
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FIG. 3. Maxwell constructions for T = 0.90 T c , T = 0.87 T c and T = 0.85 T c with the horizontal dashed line indicating the 
constant /Js/T and red triangles indicating the mixed phase boundaries at ub 1 and ns 2 . 



plotted in Fig.[2]for the highest and lowest temperatures. 
We see that all of them are more than 3cr above zero. 
This result is better than the previous ones based on the 
winding number expansion method |16l 118] , presumably 
due to the adoption of the exact projection of the deter- 
minant [19]. Thus, we believe that the sign fluctuations 
are not a problem for this study. 

We would like to point out the difference between the 
phase diagram in the grand canonical ensemble and the 
one in the canonical ensemble. We plot the expected 
canonical ensemble phase diagram in Fig. [4] in contrast 
to that in the grand canonical ensemble in Fig. [T] The 
first order phase transition line in the grand canonical 
T — p diagram becomes a phase coexistence region in the 
T - p diagram of the canonical ensemble, which has two 
boundaries that separate it from the pure phases. The 
two boundaries will eventually meet at one point. This 
point is the critical point at nonzero baryon chemical 
potential. 

Once one enters the coexistence region in a finite vol- 
ume, the contribution from the surface tension causes the 
appearance of a "double-well" in the effective free energy 
whose derivative with respect to density leads to an S- 
shaped behavior in the chemical potential versus baryon 
number plot [52]. However, in the thermodynamic limit, 
the surface tension contribution goes away since it is a 
surface term while the free energy scales with the vol- 
ume; the chemical potential will then stay constant in 
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FIG. 4. Schematic plot illustrating the scanning we use to lo- 
cate the boundaries of the mixed phase for QCD with Nf = 3. 
The infinite volume expectation for chemical potential as a 
function of density is shown in the inset. 



the coexistence phase region. The behavior of the baryon 
chemical potential in the thermodynamic limit is shown 
as an inset in Fig. [4] p\ and P2 mark the lower and upper 
boundaries of the coexistence phase at a given tempera- 
ture below T c . 

Our results for the baryon chemical potential are pre- 
sented in Fig. [5] for four different temperatures below 
T c . Statistical errors are estimated from the jackknife 
method, ft is clear that the chemical potential exhibits 
an "S-shaped" wiggle for ub between 6 and 14. To iden- 
tify the boundaries of the mixed-phase region and the 
coexistence baryon chemical potential, we rely on the 
Maxwell construction: the coexistence chemical poten- 
tial pb is the one that produces equal areas between the 
curve of the chemical potential ps as a function of Hb 
and the constant ps line which intersects with ps at «b 1 
and Ub 2 - This procedure was used in studies with stag- 
gered fermions |14[ I15j and Wilson-clover fermions |16] 
in this context for the Nj = 4 case. 

We carried out the Maxwell constructions for the three 
temperatures at 0.85 T c , 0.87 T c and 0.90 T c . We could 
not do it for the 0.93 T c case, as the wiggle there, if 
present, is not statistically significant. The results are 
presented in Fig. [3] Having determined ns ± and tib 2 for 
three temperatures, we plot the boundaries of the coex- 
istence region and perform an extrapolation in ub and T 
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FIG. 5. Phase scan for temperatures 0.85 T c , 0.87 T c , 0.90 T c 
and 0.93 T c . 
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FIG. 6. Phase boundaries in the canonical ensemble. 
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FIG. 7. Phase transition line in the T, /i plane. 



to locate the intersection of the two boundaries. To de- 
termine the crossing point, we perform a simultaneous fit 
of the boundary lines using a even polynomial in baryon 
density. We use an even polynomial since Zc is an even 
function of k. The phase boundaries and their extrap- 
olations are plotted in Fig. [6| We find the intersection 
point at T E {n%)/T c = 0.927(5) and nf = 5.7(3). 

Using the coexistence chemical potential, one can map 
out the phase diagram in the grand canonical ensemble 
as shown in Fig. [7] Note that, the region of coexistence 
phase becomes a curved transition line separating two the 
phases as we expected. In this way, we locate the critical 
point in the grand canonical ensemble at critical temper- 
ature Te/T c = 0.927(5) and baryon chemical potential 
fig/T c = 2.60(8). Using the lattice spacing a ss 0.3 fm in 



our simulation, we convert its location in physical units 
to be T E « 157 MeV and /if « 441 MeV. 

In conclusion, we have applied a canonical ensem- 
ble algorithm previously tested on the Nf — 4 to the 
more relevant Nf — 3 case and located the first order 
phase transition as signaled by the S-shapc structure 
in the /i — ng plane for several temperatures below T c . 
The Maxwell construction was employed to identify the 
boundaries of the coexistence phase and we extrapolated 
them to locate the critical point at Te — 0.925(5) T c 
and /if = 2.60(8) T c . We should point out that the 
present work is carried out on a relatively small volume 
with spatial extent of ~ 1.8 fm and for three degener- 
ate quark flavors with their masses similar to that of the 
strange quark. Quark mass for this system acts like the 
magnetic field for spin systems which weakens the phase 
transition. Since the /x = finite temperature transition 
is first order for massless quarks |23j and the present 
critical point is at a relatively large /if for quark masses 
around the strange, one expects that the critical point for 
the more realistic 2 + 1 flavor case with light u/d quarks 
to be somewhere in between. This expectation is based 
on the assumption that there is a critical surface which 
grows continuously from the critical line at /i = into fi- 
nite \x. The critical line is the one that separates the first 
order phase transition region at some finite temperature 
with small quark masses and the crossover region with 
intermediate masses (including the physical ones) on the 
fi = plane of the Columbia plot. This assumption is 
challenged by recent studies of the critical surface near 
the critical line (See Fig.l in references [51 [21]) which 
suggest that the first order region shrinks with increas- 
ing chemical potential and, therefore, there might not be 
a critical point for physical quark masses. To address 
this issue, future simulations will study the quark mass 
dependence of the critical point and the existence of the 
critical point needs to be checked on lattices with higher 
cutoffs and larger volumes. 
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